** Replication Code for main analyses for "Triangulating the relationship between education and attitudes towards immigration"
** Qinya Feng
** Journal of Experimental Political Science

/*
NOTE:
This do file performs main analysis, produces figures shown in the main text and corresponding regression results reported in the Appendix. It uses EduIA.dta generated from Data_EduIA.do.
It requires installing the STATA package reghdfe (version 5.7.3) to run fixed effects models. For further details, see http://scorreia.com/software/reghdfe/ 
The code contains three sections in the order of the studies presented in the main text and Appendix.

Sections:
1. Study 1: Discordant MZ twin analysis (Figure 2; Table A.7 - A.8)
2. Study 2: Education reform (Figure 3; Table A.13 - A.16)
3. Study 3: DZ twins & EA PGI (Figure 4; Table A.17 - A.20)
*/

clear 
cd "H:\Dofiles\Edu_IA"

set scheme s1mono
capture log close 
log using "analysis_logfile.txt", replace


use EduIA.dta, clear



**** Create variable lists for later use
global dvlist "refugee_r labor_immigration support_culture language_r"
global outreglist "treatment education_year tertiary_edu PGI_EA LOC"  
global kommunvariables "avg_RostBer avg_Valdelt avg_Andel_m avg_Andel_fp avg_Andel_c avg_Andel_s avg_Andel_v sd_RostBer sd_Valdelt sd_Andel_m sd_Andel_fp sd_Andel_c sd_Andel_s sd_Andel_v"





**********************************************
**** Study 1: Discordant MZ twin analysis ****
**********************************************
/*
The following code produce 
1) Figure 2 in the main text, 
2) Table A.7 and A.8 in the Appendix for Study 1
*/


**** The function MZ_models 
//1) runs all models in Study 1, 
//2) stores the estimates for making Figure 2, 
//3) and exports results in table formats shown in Appendix Table A.7.

** Generate an indicator for complete MZ twin pairs (necessary for keeping sample constant across model specifications)
reghdfe im_avg education_year if Zyg == 1, absorb(old_BirthKommun LopNrParID) vce(cluster LopNrParID) 
gen sampleMZ=.
replace sampleMZ=1 if e(sample)
save EduIA.dta, replace

** The function
capture program drop MZ_models
program define MZ_models

foreach var in education_year tertiary_edu{

** Models shown in the main text  
* Naive models with MZ twins
reghdfe `1' `var' i.age_2009##i.female i.edu_orientation if sampleMZ == 1, absorb(old_BirthKommun) cluster(LopNrParID)
estimate store MZ1_`1'_`var'
outreg2 using "Results\TableA7_`3'.tex", append ctitle("`2'") label keep($outreglist i.edu_orientation) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex##Age, YES, Orientation, YES, Twin FE, NO)
  
* Twin-pair fixed effect models with MZ twins; since MZ twins are the same in terms of sex and age, there is no need to control these variables here
reghdfe `1' `var' i.edu_orientation if sampleMZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
estimate store MZ2_`1'_`var'
outreg2 using "Results\TableA7_`3'.tex", append ctitle("`2'") label keep($outreglist i.edu_orientation) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex##Age, /, Orientation, YES, Twin FE, YES)

** Models only shown in the Appendix, which exclude orientation as control variables
* Naive models with MZ twins
reghdfe `1' `var' i.age_2009##i.female if sampleMZ == 1, absorb(old_BirthKommun) cluster(LopNrParID)
outreg2 using "Results\TableA7_`3'.tex", append ctitle("`2'") label keep($outreglist i.edu_orientation) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex##Age, YES, Orientation, NO, Twin FE, NO)
  
* Twin-pair fixed effect models with MZ twins; since MZ twins are the same in terms of sex and age, there is no need to control these variables here
reghdfe `1' `var' if sampleMZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA7_`3'.tex", append ctitle("`2'") label keep($outreglist i.edu_orientation) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex##Age, /, Orientation, NO, Twin FE, YES)

}

end
capture erase "Results\TableA7_imavg.tex"
capture erase "Results\TableA7_imirt.tex"


**** Table A.7 (estimates generated for Figure 2)
MZ_models im_avg "Avg index" "imavg"
MZ_models im_irt "IRT index" "imirt"
capture erase "Results\TableA7_imavg.txt"
capture erase "Results\TableA7_imirt.txt"


**** Figure 2 in the main text
capture graph drop Eduyear
capture graph drop Edutertiary
set graphics off
coefplot (MZ1_im_avg_education_year, mcolor(gs8) ciopts(color(gs8))) (MZ2_im_avg_education_year, mcolor(gs1) ciopts(color(gs1%80))), bylabel("Avg index") /// 
       || MZ1_im_irt_education_year MZ2_im_irt_education_year, bylabel("IRT index") /// 
	   ||, drop(_cons) keep(education_year) xscale(range(-0.1 0.7)) byopts(compact rows(1)) ///  
	   ylabel(, labsize(large) angle(90)) ///
	   xlabel(0.00 0.30 0.60, format(%9.2f) labsize(large))  xline(0, lpattern(dot)) ///
	   mlabel mlabsize(medlarge) format(%9.3f) mlabposition(12) legend(size(medium)) ///
	   plotlabels("Between family" "Within family") name(Eduyear)
	   
coefplot (MZ1_im_avg_tertiary_edu, mcolor(gs8) ciopts(color(gs8))) (MZ2_im_avg_tertiary_edu, mcolor(gs1) ciopts(color(gs1%80))), bylabel("Avg index") /// 
	   || MZ1_im_irt_tertiary_edu MZ2_im_irt_tertiary_edu, bylabel("IRT index") /// 
	   ||, drop(_cons) keep(tertiary_edu) byopts(compact rows(1))  ///  
	   ylabel(, labsize(large) angle(90)) ///
	   xlabel(0.00 0.30 0.60, format(%9.2f) labsize(large)) xline(0, lpattern(dot)) ///
	   mlabel mlabsize(medlarge) format(%9.3f) mlabposition(12) legend(size(medium)) ///
	   plotlabels("Between family" "Within family") name(Edutertiary)

set graphics on
graph combine Eduyear Edutertiary, ///
title("Within-family estimates for education on immigration attitudes (MZ twins)", size(medium))
graph export "Results\Figure2.pdf", as(pdf) replace


**** Table A.8: results using different items as the dependent variables
capture erase "Results\TableA8.tex"
foreach var in $dvlist{

reghdfe `var' education_year i.edu_orientation if sampleMZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA8.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex##Age, /, Orientation, YES, Twin FE, YES)

reghdfe `var' tertiary_edu i.edu_orientation if sampleMZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA8.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex##Age, /, Orientation, YES, Twin FE, YES)

}
capture erase "Results\TableA8.txt"





***********************************
**** Study 2: Education reform ****
***********************************
set matsize 2000

/*The following code produce
1) Figure 3. in the main text.
2) Table A.13, A.14, A.15, A.16 in the Appendix.
*/


**** The function reform_models 
//1) estimates the reform effect using different window widths for sample restriction
//2) stores estimates for making Figure 3.
//3) export full results for Table A.13 and A.14 in the Appendix

capture program drop reform_models
program define reform_models

local varlabel: variable label `1'

forvalues i = 10(-1)6{ //reform window widths = 10, 9, 8, 7, 6 years 

capture drop tempsample`i'

* Generate a temporary sample indicator for each window width to make sample size constant across model specifications (reghdfe iterately drop singleton observations in multiple FE groupings)
reghdfe `1' treatment i.female if sample`i' == 1 & education_year!=., absorb(BirthYear firstcohort60 reformKommun60) vce(cluster reformKommun60) 
gen tempsample`i' =.
replace tempsample`i' = 1 if e(sample)

* Export results based on the sample with 10-year window width; and generate sampleReform indicator used in the next steps 
if `i' == 10{

reghdfe `1' treatment i.female if tempsample`i' == 1, absorb(BirthYear) vce(cluster reformKommun60) //birth cohort FE included 
outreg2 using "Results\\`2'.tex", replace title(Reform and "`varlabel'") ctitle("`i'yr") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, NO, Municipality FE, NO)

reghdfe `1' treatment i.female if tempsample`i' == 1, absorb(BirthYear firstcohort60) vce(cluster reformKommun60) //reform cohort FE included
outreg2 using "Results\\`2'.tex", append ctitle("`i'yr") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, YES, Municipality FE, NO)

reghdfe `1' treatment i.female if tempsample`i' == 1, absorb(BirthYear reformKommun60) vce(cluster reformKommun60) //residential municipality FE included
outreg2 using "Results\\`2'.tex", append ctitle("`i'yr") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, NO, Municipality FE, YES)
estimate store reform1_`1'`i'

if "`1'" == "im_avg"{

* Generate a sample indicator with window width at 10-year for use in the next steps
gen sampleReform =.
replace sampleReform = 1 if tempsample10 == 1

save EduIA.dta, replace

}

}

* Export results based on 9, 8, 7, 6-year window width for sample restriction
else{

reghdfe `1' treatment i.female if tempsample`i' == 1, absorb(BirthYear) vce(cluster reformKommun60) //birth cohort FE included
outreg2 using "Results\\`2'.tex", append ctitle("`i'yr") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, NO, Municipality FE, NO)

reghdfe `1' treatment i.female if tempsample`i' == 1, absorb(BirthYear firstcohort60) vce(cluster reformKommun60) //reform cohort FE included
outreg2 using "Results\\`2'.tex", append ctitle("`i'yr") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, YES, Municipality FE, NO)

reghdfe `1' treatment i.female if tempsample`i' == 1, absorb(BirthYear reformKommun60) vce(cluster reformKommun60) //residential municipality FE included
outreg2 using "Results\\`2'.tex", append ctitle("`i'yr") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, NO, Municipality FE, YES)

}

capture drop tempsample`i'

}

capture erase "Results\\`2'.txt"

end

**** Table A.14 (estimates generated for Figure 3)
reform_models im_avg "TableA14_imavg"
reform_models im_irt "TableA14_imirt"


**** Figure 3. in the main text (sample restriction used the 10-year window width)	
coefplot (reform1_im_avg10, mcolor(gs8) ciopts(color(gs8))) (reform1_im_irt10, mcolor(gs1) ciopts(color(gs1%80))), bylabel(Reform exposure) /// 
	   ||, drop(_cons) keep(treatment) bycoefs byopts(xrescale) ///
	   ylabel(, labsize(medium) angle(90)) ///
	   xlabel(, format(%9.2f) labsize(medium)) xline(0, lc(gs3) lpattern(dot)) ///
	   mlabel mlabsize(medsmall) format(%9.3f) mlabposition(12) legend(size(medium)) ///
	   plotlabels("Avg index" "IRT index") title("Compulsory education reform and immigration attitudes", size(medium))
graph export "Results\Figure3.pdf", as(pdf) replace


**** Table A.13
reform_models education_year "TableA13"


**** Table A.15: add municipality specific pre-reform variables in the models
//These models add controls on the mean and sd of pre-reform municipality characteristics during 1940-1948; municipality fixed effects are not included. 
capture program drop reform_models2
program define reform_models2

local varlabel: variable label `1'

* Generate a temporary sample indicator (as the sample size of models with education as the dependent variable is different)
reghdfe `1' treatment i.female if sample10 == 1 & education_year!=., absorb(BirthYear firstcohort60 reformKommun60) vce(cluster reformKommun60) 
gen tempsample`1' =.
replace tempsample`1' = 1 if e(sample)

* Export results
reghdfe `1' treatment i.female `2' if tempsample`1' == 1, absorb(BirthYear) vce(cluster reformKommun60) 
outreg2 using "Results\\`3'.tex", append ctitle("`varlabel'") label keep($outreglist `2') dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, NO, Municipality FE, NO)

reghdfe `1' treatment i.female `2' if tempsample`1' == 1, absorb(BirthYear firstcohort60) vce(cluster reformKommun60)
outreg2 using "Results\\`3'.tex", append ctitle("`varlabel'") label keep($outreglist `2') dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, YES, Municipality FE, NO)

capture drop tempsample`1'
end
capture erase "Results\TableA15.tex"
reform_models2 im_avg "$kommunvariables" "TableA15"
reform_models2 im_irt "$kommunvariables" "TableA15"
//Append models with education as the dependent variable
reform_models2 education_year "$kommunvariables" "TableA15"
capture erase "Results\TableA15.txt"


**** Table A.16: results by attitude items (sample restriction used the 10-year window width)
capture erase "Results\TableA16.tex"
foreach var in $dvlist{

reghdfe `var' treatment i.female if sampleReform == 1, absorb(BirthYear) vce(cluster reformKommun60) 
outreg2 using "Results\TableA16.tex", append ctitle("`var'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, NO, Municipality FE, NO)

reghdfe `var' treatment i.female if sampleReform == 1, absorb(BirthYear firstcohort60) vce(cluster reformKommun60) 
outreg2 using "Results\TableA16.tex", append ctitle("`var'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, YES, Municipality FE, NO)

reghdfe `var' treatment i.female if sampleReform == 1, absorb(BirthYear reformKommun60) vce(cluster reformKommun60)
outreg2 using "Results\TableA16.tex", append ctitle("`var'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Birth Cohort FE, YES, Reform Cohort FE, NO, Municipality FE, YES)

}
capture erase "Results\TableA16.txt"



************************************
**** Study 3: DZ twins & EA PGI ****
************************************

/*The following code produces
1) Figure 4 in the main text.
2) Table A.17, A.18, A.19, A.20 in the Appendix.
*/


**** The function DZ_models 
//1) runs all within-family models in Study 3, 
//2) stores the estimates for making Figure 4, 
//3) and exports full results in table formats shown in Appendix Table A.17.

** Generate an indicator for complete DZ twin pairs 
qui: reghdfe im_avg PGI_EA education_year if Zyg == 2, absorb(LopNrParID) vce(cluster LopNrParID)
gen sampleDZ=.
replace sampleDZ=1 if e(sample)
save EduIA.dta, replace 

capture program drop DZ_models
program define DZ_models

* Only education 
eststo: reghdfe `1' i.female education_year if sampleDZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
estimate store DZ1_`1'
outreg2 using "Results\\`3'.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

* Only EA PGI
eststo: reghdfe `1' i.female PGI_EA if sampleDZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
estimate store DZ2_`1'
outreg2 using "Results\\`3'.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

* Both education and EA PGI included
eststo: reghdfe `1' i.female PGI_EA education_year if sampleDZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
estimate store DZ3_`1'
outreg2 using "Results\\`3'.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

end

**** Table A.17 (estimates generated for Figure 4)
capture erase "Results\TableA17.tex"
DZ_models im_avg "Avg index" "TableA17"
DZ_models im_irt "IRT index" "TableA17"
//Appending a model with education as the dependent variable
reghdfe education_year i.female PGI_EA if sampleDZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA17.tex", append ctitle("Education") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)
capture erase "Results\TableA17.txt"


**** Figure 4 in the main text
capture graph drop Edu
capture graph drop EA
set graphics off

coefplot (DZ1_im_avg, mcolor(gs8) ciopts(color(gs8))) (DZ1_im_irt, mcolor(gs1) ciopts(color(gs1%80))), keep(education_year) xscale(range(0 .2)) bylabel(Education only) ///
       || DZ3_im_avg DZ3_im_irt, keep(education_year) xscale(range(0 .2)) bylabel(Education + EA PGI) ///
	   ||, drop(_cons) byopts(compact cols(1)) ///
	   subtitle(, size(large)) ///
	   ylabel(, labsize(large) angle(90)) ///
	   xlabel(0.00 0.10 0.20, format(%9.2f) labsize(large)) xline(0, lpattern(dot) lc(black)) ///
	   mlabel mlabsize(medlarge) format(%9.3f) mlabposition(12) legend(size(medium)) ///
	   plotlabels("Avg index" "IRT index") name(Edu)

	   
coefplot (DZ2_im_avg, mcolor(gs8) ciopts(color(gs8))) (DZ2_im_irt, mcolor(gs1) ciopts(color(gs1%80))), keep(PGI_EA) xscale(range(0 .2)) bylabel(EA PGI only) /// 
	   || DZ3_im_avg DZ3_im_irt, keep(PGI_EA) xscale(range(0 .2)) bylabel(Education + EA PGI) ///
	   ||, drop(_cons) byopts(compact cols(1)) ///
	   subtitle(, size(large)) ///
	   ylabel(, labsize(large)) ///
	   xlabel(0.00 0.10 0.20, format(%9.2f) labsize(large)) xline(0, lpattern(dot) lc(black)) ///
	   mlabel mlabsize(medlarge) format(%9.3f) mlabposition(12) legend(size(medium)) ///
	   plotlabels("Avg index" "IRT index") name(EA)
	   
set graphics on
graph combine EA Edu, ///
title("Within-family estimates for EA PGI and education on immigration attitudes", size(medium))
graph export "Results\Figure4.pdf", as(pdf) replace


** Compare coefficient changes between models
//The postestimation command, suest (seemingly unrelated estimation), is used to compare coefficients across model. The suest command can adjust VCE with clustering, with stored estimates from the regress command. suest is not directly applicable to stored results from the reghdfe command. 
//Here I used the regress command to run the main models again, and use suest (with vce(cluster LopNrParID) to adjust for clustered SEs) to test whether coefficient differences are significant. 
//See STATA manual on suest usage, and Trenton et al. (2019) (DOI: 10.1177/0081175019852763) for Seemingly unrelated estimation (suest)  
capture program drop DZ_models_test
program define DZ_models_test

qui: eststo: reg `1' i.female education_year if sampleDZ == 1, absorb(LopNrParID) 
estimate store DZ1test_`1'

qui: eststo: reg `1' i.female PGI_EA if sampleDZ == 1, absorb(LopNrParID) 
estimate store DZ2test_`1'

qui: eststo: reg `1' i.female PGI_EA education_year if sampleDZ == 1, absorb(LopNrParID)
estimate store DZ3test_`1'

end
DZ_models_test im_avg
DZ_models_test im_irt

//The coefficient changes for education when including EA PGI are not significant
suest DZ1test_im_avg DZ3test_im_avg, vce(cluster LopNrParID)
lincom [DZ1test_im_avg_mean]:education_year - [DZ3test_im_avg_mean]: education_year //p=0.321 
suest DZ1test_im_irt DZ3test_im_irt, vce(cluster LopNrParID)
lincom [DZ1test_im_irt_mean]:education_year - [DZ3test_im_irt_mean]: education_year //p=0.296

//The coefficient changes for EA PGI when including education (as a mediator) are significant at the 10% level
suest DZ2test_im_avg DZ3test_im_avg, vce(cluster LopNrParID)
lincom [DZ2test_im_avg_mean]:PGI_EA - [DZ3test_im_avg_mean]: PGI_EA //p=0.088
suest DZ2test_im_irt DZ3test_im_irt, vce(cluster LopNrParID)
lincom [DZ2test_im_irt_mean]:PGI_EA - [DZ3test_im_irt_mean]: PGI_EA //p=0.091


**** Table A.18: results by immigration attitude items
capture erase "Results\TableA18.tex"
DZ_models refugee_r "refugee" "TableA18"
DZ_models labor_immigration "labor" "TableA18"
DZ_models support_culture "culture" "TableA18"
DZ_models language_r "language" "TableA18"
capture erase "Results\TableA18.txt"


**** Table A.19: between-family model
capture program drop DZ_models_between
program define DZ_models_between

eststo: reghdfe `1' i.female##i.BirthYear education_year PC1-PC20 if sampleDZ == 1, absorb(batch) vce(cluster LopNrParID)
estimate store DZ1_`1'
outreg2 using "Results\TableA19.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Twin FE, NO)

eststo: reghdfe `1' i.female##i.BirthYear PGI_EA PC1-PC20 if sampleDZ == 1, absorb(batch) vce(cluster LopNrParID)
estimate store DZ2_`1'
outreg2 using "Results\TableA19.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Twin FE, NO)

eststo: reghdfe `1' i.female##i.BirthYear PGI_EA education_year PC1-PC20 if sampleDZ == 1, absorb(batch) vce(cluster LopNrParID)
estimate store DZ3_`1'
outreg2 using "Results\TableA19.tex", append ctitle("`2'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Twin FE, NO)

end
capture erase "Results\TableA19.tex"
DZ_models_between im_avg "Avg index"
DZ_models_between im_irt "IRT index"
//Appending a model where education is the dependent variable
reghdfe education_year i.female##i.BirthYear PGI_EA PC1-PC20 if sampleDZ == 1, absorb(batch) vce(cluster LopNrParID)
outreg2 using "Results\TableA19.tex", append ctitle("Education") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Controls, YES, Twin FE, NO)
capture erase "Results\TableA19.txt"


**** Table A.20: test for Locus of control as a potential mediator
* Generate a sample indicator for complete DZ twin pairs
capture drop sampleDZ_mediator
eststo: reghdfe LOC i.female PGI_EA if sampleDZ == 1, absorb(LopNrParID) vce(cluster LopNrParID)
gen sampleDZ_mediator=. 
replace sampleDZ_mediator = 1 if e(sample)
outreg2 using "Results\TableA20.tex", replace ctitle("LOC") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

* Export model results
capture program drop DZ_models_mediator
program define DZ_models_mediator

eststo: reghdfe `1' i.female `2' if sampleDZ_mediator == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA20.tex", append ctitle("`3'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

eststo: reghdfe `1' i.female PGI_EA if sampleDZ_mediator == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA20.tex", append ctitle("`3'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

eststo: reghdfe `1' i.female PGI_EA education_year if sampleDZ_mediator == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA20.tex", append ctitle("`3'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

eststo: reghdfe `1' i.female PGI_EA `2' if sampleDZ_mediator == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA20.tex", append ctitle("`3'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

eststo: reghdfe `1' i.female PGI_EA education_year `2' if sampleDZ_mediator == 1, absorb(LopNrParID) vce(cluster LopNrParID)
outreg2 using "Results\TableA20.tex", append ctitle("`3'") label keep($outreglist) dec(3) alpha(0.001, 0.01, 0.05) addtext(Sex, YES, Twin FE, YES)

end
capture erase "Results\TableA20.tex"
DZ_models_mediator im_avg LOC "Avg index"
DZ_models_mediator im_irt LOC "IRT index"
capture erase "Results\TableA20.txt"


** Beta coefficients for all main models shown in the figures  
esttab reform1_im_avg10 reform1_im_irt10 ///
MZ2_im_avg_education_year MZ2_im_irt_education_year MZ2_im_avg_tertiary_edu MZ2_im_irt_tertiary_edu ///
DZ1_im_avg DZ1_im_irt DZ2_im_avg DZ2_im_irt DZ3_im_avg DZ3_im_irt, label beta 


clear
		
